# 铁路货运量的时序建模

# 模型评估指标
evaluate <- function(err, t, k=1, true=NA){
  mse = mean(err ** 2)
  # mse = var(err)
  rmse = sqrt(mse)
  mae = mean(abs(err))
  mape = mean(abs(err / true))
  
  aic = t * log(mse) + 2 * (k + 2)
  aicc = aic + 2 * (k + 2) * (k + 3) / (t - k - 3)
  bic = t * log(mse) + (k + 2) * log(t)
  
  print(paste('mse', mse, sep = ':'))
  print(paste('rmse', rmse, sep = ':'))
  print(paste('mae', mae, sep = ':'))
  print(paste('mape', mape, sep = ':'))
  
  print(paste('aic', aic, sep = ':'))
  print(paste('aicc', aicc, sep = ':'))
  print(paste('bic', bic, sep = ':'))
}


# 设置工作路径
setwd("F:/code/time_series/data")

# load data & 生成时序数据 ##################################################
library(readxl)

# data <- read.xlsx(file = "运输量.xlsx",sheetIndex = 2, encoding = "UTF-8")
sheet_data <- read_excel("运输量.xlsx", sheet=2)

# 查看数据类型
mode(sheet_data)


# list to frame
df = as.data.frame(sheet_data)
# 获取指定列
df_new = df[, c(1, 3)]

# 列名重置
colnames(df_new)[2] <- "freight_volume"

# 日期格式转换
df_new$month <- as.Date(paste0(df_new$month, 1), "%Y年%m月%d")
# 排序
df_new <- df_new[order(df_new$month),]

# frame to ts 格式转换
series <- ts(df_new[[2]], 
               start = c(format(start_date, "%Y"), format(start_date, "%m")), 
               end = c(format(end_date, "%Y"), format(end_date, "%m")),
               frequency = 12)
print(series)

###  切片，分训练集和测试集
series_train = window(series, start = 2015, end = 2019+6/12)
series_test = window(series, start = 2019+8/12)
print(series_train)
print(series_test)

#  数据探索 #################################################################
# 画图
library(fpp2)

autoplot(series_train) + 
  xlab("日期") + 
  ylab("铁路货运量当期值(万吨)") + 
  ggtitle("2015-2019年全国铁路货运量(万吨)") +
  theme(text = element_text(family = "STHeiti")) + 
  theme(plot.title = element_text(hjust = 0.5))

#  建模 & 评估  #############################################################
# 简单线性回归 ####################
y <- series_train
x <- seq(from=1, to=length(series_train), by=1)
line_model_1 <- lm(y ~ x)

line_train_pred = predict(line_model_1, data.frame(x=seq(from=1, to=length(series_train), by=1)))
line_test_pred = predict(line_model_1, data.frame(x=seq(from=length(series_train) + 1, to=length(series_train) + length(series_test), by=1)))

series_train_pred  <- ts(train_pred, 
                        start = c(format(start_date, "%Y"), format(start_date, "%m")),
                        end = c(2019, 8),
                        frequency = 12)

line_series_test_pred  <- ts(line_test_pred, 
                        start = c(2019, 9), 
                        end = c(format(end_date, "%Y"), format(end_date, "%m")),
                        frequency = 12)
# print(train_pred)

autoplot(series) + 
  autolayer(series_train_pred, series = "线性回归训练") +
  autolayer(series_test_pred, series = "线性回归预测") +
  xlab("日期") + 
  ylab("铁路货运量当期值(万吨) & 线性回归") + 
  ggtitle("铁路货运量简单线性回归") +
  theme(text = element_text(family = "STHeiti")) + 
  theme(plot.title = element_text(hjust = 0.5))

# 检查残差
print(checkresiduals(line_model_1))


# 模型训练评估：
line_train_err = line_series_test_pred - series_train
evaluate(err = line_err, t=length(series_train))

# 模型test评估：
line_test_err = line_series_test_pred - series_test
evaluate(err = line_test_err, t=length(series_test), true = series_test)



### 趋势性，季节性，周期性
#  季节强度

# 差分检验/平稳性

# x11 ####################
library(seasonal)
x11_series_train <- seas(series_train, x11="")
print(x11_series_train)
# 提取季节项
print(seasonal(x11_series_train))
# 提取趋势-周期项
print(trendcycle(x11_series_train))
# 提取残差项
print(remainder(x11_series_train))
# 季节调整后的时间序列
print(seasadj(x11_series_train))

autoplot(series_train) + 
  autolayer(predict(x11_series_train), series = "x11分解预测") +
  xlab("日期") + 
  ylab("铁路货运量当期值(万吨)") + 
  ggtitle("全国铁路货运量 x11分解预测") +
  theme(text = element_text(family = "STHeiti")) + 
  theme(plot.title = element_text(hjust = 0.5))



# 模型评估
summary(x11_series_train)

# 模型训练评估：


# seats ####################
sea_series_train <- seas(series_train)

autoplot(series_train) + 
  autolayer(predict(sea_series_train), series = "SEATS预测") +
  xlab("日期") + 
  ylab("铁路货运量当期值(万吨)") + 
  ggtitle("全国铁路货运量SEATS预测") +
  theme(text = element_text(family = "STHeiti")) + 
  theme(plot.title = element_text(hjust = 0.5))

# 模型评估
summary(sea_series_train)

# test评估


# stl ####################
stl_series_train <- stl(series_train, 
                    t.window = 13, 
                    # s.window = "periodic", 
                    s.window = 13, 
                    robust = TRUE)

autoplot(series_train) + 
  autolayer(predict(stl_series_train), series = "STL预测") +
  xlab("日期") + 
  ylab("铁路货运量当期值(万吨)") + 
  ggtitle("全国铁路货运量STL预测") +
  theme(text = element_text(family = "STHeiti")) + 
  theme(plot.title = element_text(hjust = 0.5))

summary(stl_series_train)


## mstl
autoplot(series_train) + 
  autolayer(predict(mstl(series_train)), series = "MSTL预测") +
  xlab("日期") + 
  ylab("铁路货运量当期值(万吨)") + 
  ggtitle("全国铁路货运量MSTL预测") +
  theme(text = element_text(family = "STHeiti")) + 
  theme(plot.title = element_text(hjust = 0.5))

autoplot(series_train) + 
  autolayer(stlf(series_train), series = "MSTL预测") +
  xlab("日期") + 
  ylab("铁路货运量当期值(万吨)") + 
  ggtitle("全国铁路货运量MSTL预测") +
  theme(text = element_text(family = "STHeiti")) + 
  theme(plot.title = element_text(hjust = 0.5))

summary(mstl(series_train))

# holt-winter 阻尼季节性方法 ####################
hw_damped_model <- hw(series_train, 
                      damped = TRUE, 
                      seasonal = "multiplicative")

autoplot(series_train) + 
  autolayer(hw_damped_model, 
            #  PI=FALSE去掉区间
            series = "holt-winter阻尼季节性乘性模型") +
  xlab("日期") + 
  ylab("铁路货运量当期值(万吨)") + 
  ggtitle("全国铁路货运量holt-winter阻尼季节性乘性模型") +
  theme(text = element_text(family = "STHeiti")) + 
  theme(plot.title = element_text(hjust = 0.5))

summary(hw_damped_model)

hw_series_test_pred = forecast(hw_damped_model, h=length(series_test))
# test评估
hw_err = time(hw_series_test_pred$mean) - series_test
evaluate(err = hw_err, t=length(series_test))


# r ets
# 参数估计
fit_params <- ets(y=series_train)

autoplot(series_train) + 
  autolayer(forecast(fit_params),
            #  PI=FALSE去掉区间
            series = "etc自动参数选择") +
  xlab("日期") + 
  ylab("铁路货运量当期值(万吨)") + 
  ggtitle("全国铁路货运量etc自动参数选择") +
  theme(text = element_text(family = "STHeiti")) + 
  theme(plot.title = element_text(hjust = 0.5))

summary(fit_params)
rs_series_test_pred = forecast(fit_params, h=length(series_test))
rs_series_test_pred = time(rs_series_test_pred$mean)
# test评估
rs_err = rs_series_test_pred - series_test
evaluate(err = rs_err, t=length(series_test), true = series_test)



# 自回归模型 arma/arima
# 检查平稳性
# Ljung-Box
# Box.test(diff(series_train), lag=1, type = "Ljung-Box")
Box.test(series_train, lag=1, type = "Ljung-Box")
# kpss检验:p>0.05,不差分
library(urca)
test <- ur.kpss(series_train)
summary(test)

# arima模型
arima_model <- auto.arima(series_train)

autoplot(series_train) + 
  autolayer(forecast(arima_model),
            #  PI=FALSE去掉区间
            series = "ARIMA模型") +
  xlab("日期") + 
  ylab("铁路货运量当期值(万吨)") + 
  ggtitle("全国铁路货运量ARIMA模型") +
  theme(text = element_text(family = "STHeiti")) + 
  theme(plot.title = element_text(hjust = 0.5))

summary(arima_model)

arima_series_test_pred = forecast(arima_model, h=length(series_test))
arima_series_test_pred = time(arima_series_test_pred$mean)
# test评估
arima_err = arima_series_test_pred - series_test
evaluate(err = arima_err, t=length(series_test), true = series_test)


# 非季节性arima
no_seasonal_arima_model <- auto.arima(series_train, seasonal = FALSE)

autoplot(series_train) + 
  autolayer(forecast(no_seasonal_arima_model),
            #  PI=FALSE去掉区间
            series = "非季节ARIMA模型") +
  xlab("日期") + 
  ylab("铁路货运量当期值(万吨)") + 
  ggtitle("全国铁路货运量非季节ARIMA模型") +
  theme(text = element_text(family = "STHeiti")) + 
  theme(plot.title = element_text(hjust = 0.5))

summary(no_seasonal_arima_model)

# 遍历arima
seasonal_arima_model <- auto.arima(series_train, 
                                   stepwise = FALSE,
                                   approximation = FALSE)
# 同上
summary(seasonal_arima_model)
